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1. Introduction Mixed integer non-linear programs combine the hardness of combinatorial explo- 
sion with the non-convexity of non-linear functions. For example, the well-known optimality conditions 
developed for differentiable objective functions have no meaning when the variables are discrete. Thus, 
it is perhaps not surprising that already linear integer programming with general quadratic constraints 
is undecidable JOI- Nevertheless, when the number of variables is fixed discrete optimization problems 
often become tractable and efficient polynomial algorithms exist (e.g. QIIHIE]). It is thus natural to ask 
what is the complexity of integer non-linear optimization assuming that the number of variables is fixed? 
We study the problem 

maximize f{xi, . . . ,Xd) subject to gi{xi, . . . ,Xd) > 0, x £ Z'^. (1) 

Here /, gi are polynomials with integral coefficients. Note that all throughout the paper we assume that 
the number of variables is fixed. Here are our two contributions to the theory: 

(1) We give a classification of the computational complexity of Problem ^ according to special cases. 
Section |2 of this article presents the details, but the reader can see the classification in Table ^ New 
results are marked with letters, known results are marked with asterisks, arrows indicate implications: 



Table 1: Computational complexity of problem Q in fixed dimension. 





Type of objective function 




convex 


arbitrary 


Type of constraints 


linear polynomial 


polynomial 


Linear constraints, integer variables 


polytime (*) <^ polytime (**) 


NP-hard (a) 




it it 




Convex semialgebraic constraints, integer variables 


polytime (**) polytime (**) 


NP-hard (c) 


Arbitrary polynomial constraints, integer variables 


undecidable (b) undecidable (d) undecidable (e) 



(2) For problem (a) , that of optimizing an arbitrary integral polynomial over the lattice points of a convex 
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rational polytope with fixed number of variables, we present an algorithm to compute a sequence of upper 
and lower bounds for its optimal value. Our bounds can be used, for instance, in a branch-and-bound 
search for the optimum. We use Barvinok's algebraic encoding of the lattice points of polytopes via 
rational functions 2 . In Section |21 we prove: 

Theorem 1 . 1 Let the number of variables d be fixed. Let f{xi , . . . ,Xd) be a polynomial of maximum total 
degree D with integer coefficients, and let P be a convex rational polytope defined by linear inequalities in 
d variables. We obtain an increasing sequence of lower bounds {ifc} and a decreasing sequence of upper 
bounds {Uk} to the optimal value 

f* — maximize f{xi,X2, ■ ■ ■ ,Xd) subject to x £ P nZ''. (2) 

The bounds L^, Uk can be computed in time polynomial in k, the input size of P and f , and the maximum 
total degree D and they satisfy the inequality Uk ^ Lk < f* ■ { \/|p7TZ^ ~ !)• 

More strongly, if f is non-negative over the polytope (i.e. f{x) > for all x E P), there exists a fully 
polynomial-time approximation scheme (FPTAS) for the optimization problem ||2J). 

We conclude with examples and a brief look at the mixed integer problem. 

2. Computational Complexity Bounds All the results we present refer to the complexity model 
where the number of operations is given in terms of the input size measured in the standard binary 
encoding. The results of H. W. Lenstra Jr. imply the entry of Table marked with (*), i.e. solving 
linear integer programming problems with a fixed number of variables can be done in time polynomial in 
the size of the input. More recently, Khachiyan and Porkolab |^ have proved that in fixed dimension, the 
problem of minimizing a convex polynomial objective function over the integers, subject to polynomial 
constraints that define a convex body, can be solved in polynomial time in the encoding length of the input. 
Thus, they settled all entries marked by (**). By the natural containment exhibited by these complexity 
classes, to show the validity of the remaining entries of Table is enough to prove the following lemma: 

Lemma 2.1 (?) The problem of minimizing a degree four polynomial over the lattice points of a 
convex polygon is NP-hard (entry (a) in Table^. 

(ii) The problem of minimizing a linear form over polynomial constraints in at most 10 integer 
variables is not computable by a recursive function (entry (b) in Table^. 

Proof. (1) We use the NP-complete problem ANl on page 249 of [7|. This problem states it 
is NP-complete to decide whether, given three positive integers a, b, c, there exists a positive integer 
X < c such that x'^ is congruent with a modulo b. This problem is clearly equivalent to asking whether 
the minimum of the quartic polynomial function (a;^ — a — by)^ over the lattice points of the rectangle 

{{x,y)\ 1 < a; < c — 1, <y< if— ^} is zero or not. This settles part (1). 

(2) In 1973 Jeroslow (TUI proved a similar result without fixing the number of variables. We follow 
his idea, but resorting to a stronger lemma. More precisely our proof relies on a 1982 result which 
states that there is no recursive function that, given an integer polynomial / with nine variables, can 
determine whether / has a non-negative integer zero, in the sense that it finds an explicit zero or returns 
null otherwise. Jones paper is a strengthening of the original solution of Hilbert's tenth problem |15|. 
Now to each polynomial / in Z[xi, X2 . . . , xq] associate the ten-dimensional minimization problem 

minimize y subject to (1 — a;2, . . . , xg) = 0, (y, xi, . . . , xg) G Z^*!}. (3) 

The minimum attained by y is either zero or one depending on whether / has an integer non-negative 
solution or not. Thus part (2) is settled. □ 

3. FPTAS for Optimizing Non-Negative Polynomials over Integer Points of Polytopes 

Consider now a polynomial function / G 'L[xi,X2, • • ■ ,Xd\ of maximum total degree D and a convex 
polytope P = {x\Ax < b} where ^ is an m x c? integral matrix and b is an integral m-vector. The purpose 
of this section is to present an algorithm to generate lower and upper bounds Lk , Uk to the integer global 
optimum value of 

maximize f{xi , . . . , Xd) subject to {xi , . . . , Xd) G P H Z'*. (4) 
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We should also remark that in our algorithm the polynomial objective function / can be arbitrary (e.g. 
non-convex) . As we have seen, the optimization problem is NP-hard already for two integer variables and 
polynomials of degree four. Nevertheless we will see that, in fixed dimension and when f{x) > for all 
X G P, the algorithm gives a fully polynomial time approximation scheme or FPTAS. This means that, in 
polynomial time on the input and (1/e), one can compute a (1 — e)-approximation to the maximum. The 
algorithm we present is based on A. Barvinok's theory for encoding all the lattice points of a polyhedron 
in terms of short rational functions. See U W for all details. Lattice points are thought of as exponent 
vectors of monomials. For example, z^z^^^ encodes the lattice point (2, —11). The set of lattice points is 
represented by a Laurent polynomial: gp{z) — X^aePnZ'' From Barvinok's theory this exponentially- 
large sum of monomials gp{z) can instead be written as a polynomial-size sum of rational functions 
(assuming the dimension d is fixed) of the form: 



gp 



(5) 



iei 



n(i 



where / is a polynomial-size indexing set, and where E-i G {1,-1} and Ui,Vij £ TH^ for all i and j. For 
details see [21IS.,, 6 . 

We need a way to encode via rational functions the values of the polynomial / over all the lattice points 
in a polytope. The key idea, first introduced in Lemma 9 of 6 and generalized in 0, is that differential 
operators associated to / can be used to compute a rational function representation of X^aePnZ'' fi'^)^"'- 
The following Lemma recently appeared in ^3^ : 

Lemma 3.1 Let gp{z) be the Barvinok representation of the generating function of the lattice points of P. 
Let f be a polynomial in Z[xi, . . . , Xd] of maximum total degree D. We can compute, in time polynomial 
on D and the size of the input data, a Barvinok rational function representation gpj(z) for the generating 
function J^aePnE-^ f{a)z''. 



Proof. We give here the author's original proof the lemma for D fixed. The first proof without this 
assumption was recently given by A. Barvinok in W- 

We begin assuming f{z) = Zr, the general case will follow from it: Consider the action of the differential 
operator Zr-^ in the generating function gp{z) and on its Barvinok representation. On one hand, for 



Zr-g^ ■ gp{z) 



d 



the generating function 

QePnz<i aGPnZ'' 
On the other hand, by linearity of the operator, we have that in terms of rational functions 

f \ 



dZr 



iei 



dZr 



n(i-^"-) 



Thus it is enough to prove that the summands of the expression above can be written in terms of rational 
functions computable in polynomial time. The standard quotient rule for derivatives says that 



_d_ 

dZr 



(^)n-=i(i 



We can expand the numerator as a sum of no more than 2'^ monomials. This is a constant number 
because d, the number of variables, is assumed to be a constant. This argument completes the proof of 
our lemma when /(z) = z^. 

For the case when /(z) is a general monomial, i.e. /(z) = c • zf^ • . . . • z^'', then we can compute again 
a rational function representation of gpj{z) by repeated application of basic differential operators: 



c Zi 



d 
dzi 



zd^\ -gpiz) 



E 



c-a''z". 
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Thus we require no more than 0{D'^) repetitions of the single-variable case. 

Finally, if we deal with a polynomial / of many monomial terms, we compute and add up all such 
expressions that we get for each term of f{x) and obtain our desired short rational function representation 
for the generating function for X^aePnZ'' fi^)^"- Note that only polynomially many steps are needed 
because d is fixed and the largest number of possible monomials in / of degree s is {'^^1^^), thus for 
fixed d we will do no more than 0{D'^) repetitions of the monomial case. □ 

Now we are ready to present our algorithm to obtain bounds Uk, Lk that reach the optimum. Step 1 of 
preprocessing is necessary because we rely on the elementary fact that, for a collection S — {si, . . . , Sr} 

of non-negative real numbers, maximum{si|si G S} equals limfc^oo \J'Y^^j=i^ ■ 
Algorithm 

Input: A rational convex polytope P C M*^, a polynomial objective / € Z[a;i, . . . ,Xd\ of maximum total 
degree D. 

Output: An increasing sequence of lower bounds Lfe, and a decreasing sequence of upper bounds Uk 
reaching the maximal function value /* of / over all lattice points of P. 

Step 1. If / is known to be non- negative in all points of P, then go directly to Step 2. Else, solving 2d 
linear programs over P, we find lower and upper integer bounds for each of the variables xi, . . . , a;^. Let 
M be the maximum of the absolute values of these 2d numbers. Thus \xi\ < M for all i. Let C be the 
maximum of the absolute values of all coefficients, and r be the number of monomials of f{x). Then 

L -rCM° < f{x) < rCM° U, 

as we can bound the absolute value of each monomial of f{x) by CM^ . Replace / by f{x) — f{x) ~ L < 
U — L, a non- negative polynomial over P. Go to Steps 2, 3, etc. and return the optimal value of /. 
Trivially, if we find the optimal value of / over P we can extract the optimal value for /. 

Step 2. Via Barvinok's algorithm (see 01121 EI )j compute a short rational function expression for the 
generating function gp{z) = '^aepnz'^ From gp{z) compute the number |P n Z'^| — gp{l) of lattice 
points in P in polynomial time. 

Step 3. From the rational function representation gp{z) of the generating function ^ compute 

the rational function representation of gpjk{z) of X^aePnZ'' /'^(o^)-^" polynomial time by application 
of Lemma 13.11 We define 



Lk ygpji'W/gpjoil) and Uk --^ ygpjk{i). 

When [Uk\ — [Pfc] < 1 stop and return [L^] — [Uk\ as the optimal value. 
End of Algorithm. 



Lemma 3.2 The algorithm is correct. 



Proof. Using the fact that the arithmetic mean of a finite set of nonnegative values is at most as big 
as the maximum value, which in turn is at most as big as the sum of all values, we obtain the sequences 
of lower and upper bounds, Lk and Uk, for the maximum: 



Lk = 



\ 



E /(")' 

ctGPnZ'i 

~~\PnW\~ 



< max{/(a) : a G P n Z"^} < ^ 




fH" - Uk. 



Note that as s —>■ oo, Lk and Uk approach this maximum value monotonously (from below and above, 
respectively). Trivially, if the difference between (rounded) upper and lower bounds becomes strictly less 
than 1, we have determined the value max{/(a;) : x G P Ci Z*^} = \Lk~\- Thus the algorithm terminates 
with the correct answer. □ 



Theorem 11.11 will follow from the next lemma: 
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Lemma 3.3 Let f be a polynomial with integer coefficients and maximum total degree D. When the 
dimension d is fixed, 

(i) the bounds Lk, Uk can be computed in time polynomial in k, the input size of P and f , and the 
total degree D. The bounds satisfy the following inequality: 

c/fc-ife<r • 

[a) In addition, when f is non-negative over P (i.e. fix) > for all x G P), for fc = (1 + 
1/e) log(|P n 'Z'^D, Lk is a [1 — e)- approximation to the optimal value f* and it can be computed 
in time polynomial in the input size, the total degree D, and 1/e. Similarly, Uk gives a (1 + e)- 
approximation to f* . Moreover, with the same complexity, one can also find a feasible lattice 
point that approximates an optimal solution with similar quality. 

Proof. Part (i). From Lemma im on fixed dimension d, we can compute gpj = J^aePm.-^ fi'^)^"' 
a rational function in time polynomial in D, the total degree of /, and the input size of P. Thus, because 
f'' has total degree of Dk and the encoding length for the coefficients of f'^ is bounded by k \og{kC) 
(with C the largest coefficient in /), we can also compute gpjk — X^aePnZ'' /'^(o^)^" in time polynomial 
in fc, the total degree D, and the input size of P. Note that using residue techniques P], we can evaluate 
gp fk{l) in polynomial time. Finally observe 



Uk~Lk^ J J2 
y aePnZ'' 

= Lk (^iPnZ'il 



Part (ii). Note that if \/\P n I/\ — 1 j < e then Lk is indeed a (1 — e)-approximation because 

f* <Uk^Lk + [Uk - Lk) < Lk + f*(^ ^IPDZ^I -lj<Lk + f*e. 

Observe that (j>{e) := (1 + l/e)/(l/ log(l + e)) is an increasing function for e < 1 and lime_,o 0(e) — 1, 
thus 0(e) > 1 for < e < 1. Hence, for aU fc > log(|P n Z'']) + log(|P n Z'^D/e > log(|P n Z''])/ log(l + 

e), we have indeed ^ \/\P fl Z''! — 1^ < e. Finally, from Lemma l3.ll the calculation of Lk for fc — 

log(|P n Z'^l) + log(|P n Z'^D/e would require a number of steps polynomial in the input size and 1/e. A 
very similar argument can be written for Uk but we omit it here. 

To complete the proof of part (ii) it remains to show that not only we approximate the optimal value 
/* but we can also efficiently find a lattice point a with /(a) giving that quality approximation of /*. 
Let fc = (1 + 1/e) log(|P n Z''!), thus, by the above discussion, Lk is an (1 — e)-approximation to /*. Let 
Qo := [—M,M]'^ denote the box computed in Step 1 of the algorithm such that P C Qq. By bisecting 
Qo, we obtain two boxes Q[ and Q". By applying the algorithm separately to the polyhedra P n Q[ and 
P n Qi, we compute lower bounds LJ. and L'^ for the optimization problems restricted to Q[ and Q", 
respectively. Because is the arithmetic mean of f''{a) for a e P n Z'*, clearly 

mm{L[,Ll} < Lk < max{L',., L'^}. 

Without loss of generality, let LJ, > L'^. We now apply the bisection procedure iteratively on Q'^,. After 
dlogM bisection steps, we obtain a box Q'j, that contains a single lattice point a G P fl QJ. fl Z'^, which 
has an objective value /(a) = ij. > Lk > (1 — e)/*- C 

We remark that if we need to apply the construction of Step 1 of the algorithm because / takes negative 
values on P, then we can only obtain an (1 — e)-approximation (and (1 + e)-approximation, respectively) 
for the modified function / in polynomial time, but not the original function /. We also emphasize that, 
although our algorithm requires the computation of J2a£P fi*^) ^^r different powers of /, these numbers 
are obtained without explicitly listing all lattice points (a hard task), nor we assume any knowledge of 
the individual values /(a). We can access the power means X^agp /'(o^) indirectly via rational functions. 
Here are two small examples: 



-i) <r (^/\Pnzd\-i) . 
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Example 1, monomial optimization over a quadrilateral: The problem we consider is that of 
maximizing the value of the monomial x^y over the lattice points of the quadrilateral 

{(2;,y)|3991 < 3996 a; -4y < 3993, l/2<x< 5/2}. 

It contains only 2 lattice points. The sum of rational functions encoding the lattice points is 

x^yiooo ^ xy , xy , a;2yiooo 



(l-(xy999)-i)(l-?y-i) (1 - (1 - {1 - xy^^^) {1 - y) {1 - (xy^^^)-^) {1 - y)' 

In the first iteration Li = 4000.50 while Ui = 8001. After thirty iterations, we see L30 = 7817.279750 
while [/30 — 8000, the true optimal value. 

Example 2, nvs04 from MINLPLIB: A somewhat more complicated example, from a well-known 
library of test examples (see |http : / / www . gamswor Id . org/ minlp7| , is the problem given by 

2 




min 100 I - + - I - + ii I | + { --h 



2 



(6) 



s.t. H,i2 e [0,200] nz. 



Its optimal solution as given in MINLPLIB is ii = 1, 22 = 2 with an objective value of 0.72. Clearly, to 
apply our algorithm from page0]literally, the objective function needs to be multiplied by a factor of 100 
to obtain an integer valued polynomial. 

Using the bounds on ii and i2 we obtain an upper bound of 165 • 10^ for the objective function, which 
allows us to convert the problem into an equivalent maximization problem, where all feasible points have 
a non-negative objective value. The new optimal objective value is 164999999999.28. Expanding the new 
objective function and translating it into a differential operator yields 

4124999999947 d 172 d / d f d \ 

— 25 — " - ^"'-'ar, + a;: - air j - '"^ V'ar^ 



OZ2 J \ ozi J \ 0Z2 J \ OZl J \ OZl J \ OZl J 

The short generating function can be written as 5(21, ^2) = (^i"rj7 ~ r^ir) (1^2 1-Z2 

In this example, the number of lattice points is |P n = 40401. The first bounds are Li = 
139463892042.292155534, Ui = 28032242300500.723262442. After 30 iterations the bounds become 
L30 = 164999998845.993553019 and C/30 = 165000000475.892451381. 

4. An extension to the mixed integer case Now, we wish to discuss extensions of Theorem ll.il 

to the mixed integer scenario. If some of the variables are allowed to be continuous then we can describe 
the task as 

maximize /(x, y) subject to {(x, y)\Ax + By < h, with Xi G Z and yi G M}. 

Consider the sequence of integer polynomial optimization problems optimize f{x, — ) subject to (x,?/) G 
r?i = {{x,y)\Ax + B{^) < b, with Xi G Z and yi G Z}, where each of the subproblems is equivalent to 
optimizing a polynomial over a "semi-dilated" polytope (in some coordinate directions but not others). 
As n goes to infinity, the sequence of optimal solution values can have several limit points. Nevertheless, 
it is still possible to construct a subsequence of problems whose optimal values approximate the mixed 
integer optimum to arbitrary precision: 

Corollary 4.1 With the hypotheses of Theorera M.R we can construct a sequence of integer polynomial 
programming problems, over finer and finer grids, whose optimal values converge to the optimal value of 
the mixed integer program 

f* — maximize f{xi, X2, • • ■ , Xfi) subject to x G P, and x^ G Z for i E I C_ {1, 2, . . . , d}. 

By applying the algorithm of Theorem \l.l\ to the subproblems, we can approximate the optimum to arbi- 
trary precision. 
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When all variables are continuous, the original polytope P is actually dilated uniformly in all directions 
by a parameter n, or equivalently, the integer grid is refined. Define 

QG(nP)nZ<' 

for each power q and dilation factor n. Note that for fixed q, then we can easily see that the sum 
'^ae(nP)nZ'^ 'K^fi'^)^'^ essentially an approximation to the Riemann integral of f^; thus 

lim r{q,n)/ 71,'^ = / f'^{x)dx. 



As n, q grow, the values ^r{q, n) approximate the sequence ^ Jp fi{x)dx which converges to 

(max^:^ p/(a;)) • volume(P). This is related to recent work (see 0)^1^] and references therein) where the 
global optimum of a polynomial over a compact domain is investigated as the result of a grid refinement 
and properties of sums of squares. 
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